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FOREWORD 


The  Center  for  Air  Sea  Technology  (CAST)  research  program  in  FY96  was 
modified  to  adjust  to  new  Navy  Ocean  Modeling  and  Prediction  (NOMP) 
program  priorities,  especially  in  the  area  of  coastal  and  semi-enclosed  seas.  The 
objectives  were  to: 

•  Conduct  coastal  and  semi-enclosed  seas  ocean  modeling  basic  research, 
embedded  in  a  CAST  modularized  software  system,  with  the  emphasis 
on  model  relocatability  to  any  geographical  region; 

•  Support  the  technical  requirements  of  Navy  and  university  ocean 
modeling  efforts  by  providing  routine  day-to-day  technical  support  to 
the  scientific  staff,  and  by  designing,  developing,  and  implementing 
next  generation  technical  support  capabihty; 

•  Tailor  and  transition  applicable  advanced  technical  support  capabilities 
developed  for  the  research  community  to  the  operational  Navy;  and 

•  Strengthen  collaboration  with  academia  by  incorporating  student  and 
faculty  in  CAST  projects. 

In  accomplishing  these  objectives,  CAST  in  FY1996  supported  two  postdoctoral 
assistants.  This  technical  note  summarizes  the  FY1996  research  conducted  by 
these  individuals  and  their  advisors.  CAST  was  extremely  pleased  with  the 
research  support  provided  by  these  individuals, 
the  quality  of  the  research  conducted. 
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fflGH  RESOLUTION  COASTAL  STUDIES 
by  Dr.  Avichal  Mehra 


SCIENTIFIC  ADVISOR;  Dr.  David  Dietrich,  Center  for  Air  Sea  Technology 
OBJECTIVES: 

To  numerically  model  basic  coastal  phenomenon  using  accurate  and  efficient 
numerics  at  high  resolution,  and  interpret  them  in  terms  of  theory  and 
observations  to  achieve  a  better  understanding  of  the  complex  ocean  dynamical 
processes  involved. 

APPROACH; 

The  z-level  Mississippi  State  University  DieCAST  (Dietrich  Center  for  Air  Sea 
Technology)  Model  has  emerged  as  very  accurate  and  robust  in  application  to 
small-scale  near-coast  features  such  as  boundary  currents  and  shelf-break  fronts, 
and  their  interaction  with  deep  water  eddies.  With  the  Navy's  emphasis  on  coastal 
and  shallow  water  regions,  there  is  a  need  for  high  resolution  models  that 
realistically  and  efficiently  simulate  coastal  ocean  dynamics,  including  density 
fields  and  boundary  currents. 

The  modified  Arakawa  "a"  grid  DieCAST  Model  is  three-dimensional  in  a 
rotating  frame,  and  uses  the  hydrostatic,  Boussinesq,  and  rigid-lid  approxima¬ 
tions.  It  also  uses  flux  conserving  centered  approximations  based  on  control 
volumes  and  a  weakly  time  filtered  leapfi’og  time  integration  scheme. 

A  significant  feature  is  that  the  model  applies  interpolations  to  communicate  data 
between  collocated  "a"  and  staggered  "c"  grid  locations.  Such  mixed  grid 
approaches  mitigate  weak  points  of  the  respective  grids(Coriolis  terms  on  the  "c" 
grid;  "null  space"  pressure  problems  on  the  "a"  grid).  To  improve  grid 
convergence  and  reduce  dispersion,  higher  order  interpolations  are  needed. 
Hence,  the  two  following  modifications  were  introduced: 

1.  Reduced  Dispersion  Advection  Scheme.  The  first  modification  involves  the 
interpolations  used  for  the  advection  terms.  This  scheme  greatly  reduces 
numerical  dispersion,  and  is  called  the  reduced  dispersion  advection,  or 
"RDA",  scheme.  The  Arakawa  "a"  grid  model  uses  a  "c"  grid  non-divergent 
advection  velocity.  The  control  volume"a"  grid  advection  scheme  thus 
requires  evaluation  of  the  fluxed  quantities  at  the  staggered  "c"  grid  locations. 
Instead  of  using  the  classic  linear  advection  scheme  based  on  two-point 
averages,  used  by  the  early  versions  of  DieCAST  ocean  model,  the  RDA 
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scheme  includes  additional  infonnation  from  the  next  nearest  points  which  is 
exact  for  general  cubic  variations.  The  resulting  interpolated  field  is  then 
multiplied  by  the  non-  divergent  advection  velocity  to  get  the  flux  at  the  cell 
face.  These  fluxes  are  then  substituted  into  the  "control  volume"  integrated 
conservation  equations. 

2.  Modified  Incompressibility  Algorithm.  The  second  modification  to  the  ocean 
model  involves  the  numerical  algorithm  for  construction  of  a  non-divergent 
advection  velocity,  and  is  referred  to  as  the  modified  incompressibility 
algorithm,  or  "MIA"  scheme.  As  with  the  original  "a"  grid  scheme,  the 
non-divergent  advection  velocity  was  determined  by  interpolating  the  "a" 
grid  trial  velocity  (derived  using  a  trial  surface  pressure,  the  hydrostatic 
relation  and  the  momentum  equations)  to  the  usual  "c"  grid  advection 
locations,  and  then  "clearing  out  the  divergence"  of  the  resulting  "c"  grid 
velocity  by  appropriate  surface  pressure  adjustment  and  application. 
Originally,  this  incompressibility  information  was  communicated  back  to  the 
"a"  grid  locations  by  simply  interpolating  the  "c"  grid  advection  velocity  back 
to  the  "a"  grid  locations,  using  the  cubic  spline  approach  in  reverse.  In  order 
to  take  advantage  of  the  greatly  reduced  numerical  dispersion  of  the  RDA 
scheme  described  above,  and  thus  achieve  low  overall  numerical  dispersion, 
the  reverse  interpolation  approach  is  modified.  Instead  of  interpolating  the 
velocity,  the  changes  of  velocity,  that  result  from  "clearing  out  its 
divergence",  are  interpolated.  In  other  words,  the  interpolated  changes  are 
added  to  the  original  "a"  grid  "trial"  velocity  referred  to  above. 

Combining  this  MIA  scheme  with  the  RDA  scheme,  the  overall  numerical 
dispersion  is  greatly  reduced.  The  changes  of  velocity  needed  to  maintain 
incompressibility  are  small  compared  to  velocity  differences  at  adjacent  grid 
locations,  because  a  good  guess  for  the  "trial"  top  level  pressure  field  is  available 
each  time  step  from  previous  results,  thus  producing  a  quasi-non-divergent 
velocity  that  requires  only  a  small  adjustment  to  "clear  out  its  divergence".  These 
new  numerics  are  applicable  to  both  free  surface  models  as  well  as  general 
circulation  models.  Much  needs  to  be  done  to  understand  the  behavior  of  these 
schemes  in  detail,  but  the  results  (given  in  the  next  section)  show  that  they  have 
interesting  effects  and  are  quite  promising. 

RESULTS: 

The  above  low  dispersion  numerics  were  applied  to  two  different  studies 
involving  distinct  ocean  dynamical  processes.  The  first  was  the  semi-enclosed 
basin  of  the  Gulf  of  Mexico  which  is  one  of  the  most  thoroughly  observed  ocean 
basins.  It  has  well  known  mean  boundary  conditions,  which  control  the  internal 
Loop  Current  and  eddy  shedding  dynamics.  The  Gulf  is  large  enough  to  include 
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major  ocean  basin  dynamic  features,  yet  small  enough  to  model  with  high 
resolution  on  workstations. 

Model  parameters  included  a  lateral  resolution  of  1/12°  and  20  vertical  layers. 
Real  unfiltered  topography  was  used,  with  depths  truncated  at  3.5  km,  and 
climatological  winds  and  30-day  surface  restoring  to  Levitus  climatology  were 
applied.  Temperature-only  based  thermodynamics  were  used  for  the  GOM  study 
which  is  appropriate  when  deep  water  formation  effects  are  secondary.  Initial  and 
inflow  temperatures  were  based  on  annual  mean  Levitus  climatology.  Eastern 
open  boundary  inflow  temperatures  were  then  adjusted  to  give  surface  trapped 
geostrophic  inflow  consistent  with  the  observed  data;  the  corresponding  inflow 
velocities  were  specified  as  time-constant  inflows.  An  approximately  16  Sverdrup 
barotropic  western  boundary  current  inflow  was  specified  to  represent  inflows 
from  the  southern  Caribbean,  giving  a  total  Caribbean  inflow  of  25  Sverdrups. 
The  initial  velocity  was  zero,  but  inflows  were  switched  on  full  strength  at  the 
first  time  step. 

Lateral  heat  and  momentum  diffusivities  were  a  constant  10  m-m/sec  and  annual 
mean  Hellerman  winds  were  applied  for  surface  forcing.  Vertical  heat  and 
momentum  diffusivities  were  1  cm-cm/sec  plus  a  small  additive  variable 
numerical  contribution.  The  present  simulation  was  done  on  the  NAVOCEANO 
Cray  C-90  with  210  megaflops  computation  speed.  With  the  1/12°  lateral  resolu¬ 
tion  used  on  a  194  X  146  X  20  layers  grid,  the  model  simulated  approximately 
4  years  of  real  GOM  time  per  cpu-day. 

The  new  numerics  were  robust  in  application  to  the  propagation  of  warm  core 
Loop  Current  eddies  across  the  Gulf  of  Mexico.  This  propagation  involved  highly 
nonlinear  dynamics:  a  strong  front  was  maintained  at  the  outer  edge  of  the  warm 
core  eddy  as  it  spun  at  about  one  rotation  per  15  days  while  propagating  across 
the  Gulf  of  Mexico  over  a  period  of  about  six  months  (see  Figure  1).  Frontal 
eddies  slowly  dispersed  the  warm  core  water  as  it  spun  westward.  Compared  to 
the  earlier  model  results,  the  new  model  results  show  that  the  warm  core  eddies 
deform  more  and  lose  their  identity  more  rapidly  because  of  this  dispersion,  but 
the  dispersed  energy  appears  in  a  much  larger  number  of  eddies,  as  observed. 

The  new  modified  schemes  decrease  numerical  dispersion,  and  lead  to  tighter 
fronts  and  small  scale  frontal  eddies  (Figure  1).  These  can  mix  material  across 
fronts,  resulting  in  "physical"  (involving  resolved  velocity  components  rather 
than  interpolations)  dispersion.  Such  dispersion  can  extend  the  mixing  region  for 
different  water  masses,  leading  to  a  more  widely  distributed  eddy  field.  For  small 
scale  flow  components,  such  as  frontal  eddies,  such  "physical"  dispersion  can  of 
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Equivalent  Sea  Surface  Height  (in  cms)  on  Day  1590 
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course  have  significant  truncation  error  effects,  yet  the  very  existence  of  the 
frontal  eddies  on  tight  fronts  is  an  improvement  in  representing  nature. 

The  second  study  involved  an  idealized  lake  Kelvin  wave  problem  to  study  fast 
(non-geostrophic)  internal  modes  of  ocean  physics,  specifically  the  transient 
Poincare  and  Kelvin  wave  dynamics.  The  computational  basin  was  a  100  km 
diameter  circular  flat  bottom  lake  which  was  100  m  deep  and  patterned  after  the 
Great  Lakes.  The  domain  was  forced  by  wind  which  ramps  up  for  a  day  and  then 
shuts  off  after  two  days.  It  included  the  dominant  geophysical  dynamic  effects 
(buoyancy  and  rotation)  which  lead  to  formation  of  Kelvin  edge  waves  that  are 
challenging  to  resolve  and  maintain  against  numerical  dissipation  and  dispersion. 

A  linear  equation-of-state  was  used  with  different  grid  resolutions  of  5.00,  2.50, 
1.25  and  0.6125  kms  to  check  convergence.  The  new  reduced  dispersion  scheme 
gave  close  to  the  correct  linear  Kelvin  wave  phase  speed  (about  8  days  for 
propagation  around  the  lake)  even  at  coarse  resolutions.  The  results  had  very 
small  amount  of  dissipation  as  the  warm  front  of  the  marginally  resolved  Kelvin 
edge  wave  was  clearly  seen  for  eight  revolutions  around  the  basin  with  only 
2.5  km  resolution  (80  x  80  grid).  (An  animation  of  these  results  is  available  on 
the  web  site  at  http:\\www.cast.msstate.edu.) 

The  actual  patterns  were  found  to  be  rapidly  changing  with  time  scales  of  order 
1  to  8  days  dominating.  Initially  the  Poincare  waves  dominate  the  results  as  seen 
in  Figure  2.  But  thereafter,  the  higher  resolution  low  dispersion  numerics  cause 
the  Kelvin  wave  warm  front  to  tighten  greatly  after  day  20  (see  Figure  3).  Such 
tightening  is  consistent  with  the  fact  that  Kelvin  wave  warm  fronts  are  much 
more  commonly  observed  than  cold  fronts.  As  this  happens,  along-shore  scale  of 
the  trailing  warm  water  decreases  rapidly  until  its  along-shore  scale  roughly 
matches  its  off-shore  scale,  which  is  the  natural  Rossby  radius  of  deformation.  It 
thus  becomes  a  half-cylinder  warm  core  trapped  against  the  boundary  moving  at 
the  natural  Kelvin  wave  speed. 

FUTURE  RESEARCH; 

Plans  for  FY97  and  FY98  include  the  following: 

Task  1:  Nested  Santa  Barbara  Channel  Model:  A  new  high  resolution  Santa 
Barbara  Channel  (SBC)  model,  using  1/48  deg  (about  2  km)  resolution,  will  be 
nested  inside  an  intermediate  resolution  (1/12  deg)  California  Current  (CC) 
model  developed  in  collaboration  with  NPS.  Both  will  use  the  accurate,  reduced 
dispersion  numerics  of  the  new  DieCAST.  Time  permitting,  the  CC/SBC  system 
will  be  nested  in  the  NRL  Layered  Pacific  Model.  The  initial  coupling/nesting 
will  be  one-way-only,  from  large  to  small  model  subdomains.  Time  permitting, 
two-way  coupling  will  be  implemented  in  a  triply  nested  system. 
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Transient  Wind  Forced  Poincare  Waves 
Temperature  (in  C)  at  Depth  9.7  m  on  Day  8.00 
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Transient  Wind  Forced  Kelvin  Edge  Waves 
Temperature  (in  C)  at  Depth  9.7  m  on  Day  30.00 
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Larger  domains  will  provide  normal  boundary  velocity  (nbv)  conditions  used  by 
successively  smaller  domains  in  the  nesting  hierarchy.  The  nbv  will  be 
distributed  non-uniformly  according  to  the  nested  model  internal  dynamics,  but 
the  summed  volume  flows  across  a  coarse  resolution  interface  will  be  constrained 
to  match  the  coarse  resolution  results.  Normal  boundary  fluxes  will  be 
determined  by  the  smaller  domain  models.  A  distributed  jet  approach  may  be 
used  for  the  resulting  momentum  sources  to  coarser  grids.  Mass,  momentum, 
temperature  and  salinity  will  be  conserved  by  the  coupling  approach. 

Task  2:  Data  Assimilation:  This  work  will  be  leveraged  by  work  done  by  the 
Minerals  Management  Service  (MMS).  The  ultimate  MMS  goal  is  a  nowcast/ 
forecast  system  for  the  California  Bight  Region.  Thus,  the  work  will  include  data 
assimilation  into  an  operational  version  of  the  DieCAST  model  for  MMS.  The 
data  assimilating  DieCAST  SBC  model  will  be  validated  by  the  large  SBC  data 
base.  Promising  data  assimilation  approaches  will  be  implemented  and  evaluated, 
as  they  become  available,  in  collaboration  with  NPS,  MIT,  and  Dalhousie. 

Task  3:  General  Coastal  Model:  A  special  coastal  version  of  DieCAST  will  be 
developed  which  conforms  to  irregular  coastlines  by  varying  the  logical  grid  in  a 
way  such  that  most  coastal  land  points  will  be  eliminated.  By  shifting  indices  row 
by  row,  this  will  lead  to  a  code  with  good  vector  and  parallel  computer  efficiency 
nesting  capability,  for  an  irregular  domain  model  which  will  be  developed  later. 

Task  4:  Other  Applications:  Idealized  studies  including  the  sensitivity  of  the  basic 
phenomena  in  wakes  from  individual  islands  to  the  various  parameters  (steepness, 
Reynolds  number,  stratification,  latitude)  will  be  explored.  These  will  also 
include  island  encounters  with  eddies  advected  from  upstream  regions,  and 
topography  effects  such  as  steering  by  nearby  ridges  and  seamounts.  Also,  after 
achieving  better  understanding  of  the  basic  physics  of  island  flows,  the  flow 
dynamics  around  real  islands  in  the  Ryukyu  island  chain  will  be  studied.  The 
East  China  Sea  is  particularly  challenging  because  it  includes  a  strong  boundary 
current  along  the  continental  shelfbreak  on  its  west  side,  and  also  includes  strong 
water  mass  transformations  during  wintertime.  Accurate  simulations  requires 
high  resolution,  particularly  over  the  shelfbreak  and  island  chains  for  which  large 
computing  resources  will  be  needed.  To  reduce  the  computing  costs,  we  will 
share  the  computing  burden  by  collaborating  with  Dr.  Steve  Piacsek  who  has 
developed  a  massively  parallel  version  of  DieCAST  for  the  CM5.  An  alternative 
to  the  East  China  Sea  is  the  Eastern  Caribbean  Sea,  because  of  the  large  amount 
of  available  data,  including  a  1/16  resolution  West  Indies  data  set  prepared  by 
Drs.  T.  Townsend  and  A.  McManus  of  NRL. 
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RELOCATABLE  NUMERICAL  MODELS  OF  MARGINAL, 
SEMI-ENCLOSED,  AND  COASTAL  SEAS 
by  Dr.  Avichal  Mehra 


SCIENTIFIC  ADVISOR:  Dr.  Lakshmi  Kantha,  University  of  Colorado 

TECHNICAL  ADVISOR:  Mr.  Valentine  Anantharaj,  Center  for  Air  Sea 

Technology 

OBJECTIVE: 

To  develop  and  comprehensively  test  and  validate  increasingly  complex 
numerical  models  of  the  coastal,  marginal,  and  semi-enclosed  seas  with  the  entire 
process  integrated  into  a  relocatable  modeling  environment,  and  then  transition 
the  integrated  system  to  the  Navy's  research  and  operational  communities. 

APPROACH: 

Under  ONR/SPAWAR  sponsorship,  MSU  CAST  and  Dr.  Lakshmi  Kantha  of  the 
University  of  Colorado  have  developed  a  two-dimensional  tidal  model  embedded 
in  a  user  interactive  interface  that  allows  quick  relocation  of  the  model  anywhere 
geographically  for  which  databases  and  forcing  fields  exist.  This  model,  called 
CURReNTSS,  is  a  finite  difference,  explicit,  vertically-integrated  barotropic 
model  assimilating  tidal  component  data  from  coastal  tide  gages  and  from  any 
bottom  pressure  gages  available.  It  is  fully  non-linear  and  a  sub-  component  of 
the  NAVOCEANO/CU  3-D  operational  model.  The  model  incorporates  direct 
astronomical  forcing  and  can  also  include  surface  forcing  (winds  stress  and 
pressure  fields)  from  atmospheric  models,  if  available,  to  predict  storm  surges. 
The  open  boundary  conditions  are  obtained  from  an  earlier  global  run  at  1/5° 
resolution,  which  are  readily  accessible  from  the  GUI. 

The  GUI  also  provides  bathymetry  from  the  ET0P05  database,  which  can  be 
refined  within  the  GUI  and  setup  for  the  model.  It  is  also  possible  to  import 
bathymetry  from  other  databases  into  the  GUI.  The  model  can  be  run  with  any 
number  of  components  including  long  term  and  compound  tides.  It  employs  a 
simple  data  assimilation  procedure  by  replacing  the  model  predicted  SSH  at 
pre-determined  intervals  by  a  weighted  sum  of  the  model  prediction  and  the 
observed  SSH  from  the  tide/bottom  pressure  at  that  grid  point,  the  weights 
determined  a  priori.  The  tide  gage  data  comes  from  the  data  base  at  the 
International  Hydrographic  Organization  (1979),  supplemented  by  the  Admiralty 
charts  (1993).  This  database  is  also  available  from  the  GUI  and  relevant  tidal 
stations  can  be  easily  extracted  for  assimilation  in  the  model. 
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The  second  real-time  case  study  was  carried  out  in  April  1996,  in  response  to 
political  turmoil  in  the  West  AJErican  nation  of  Liberia.  Once  again,  the  study 
was  carried  out  at  two  different  resolution  in  a  nested  fashion.  The  physical 
domain  for  the  coarse  grid  stretched  from  the  equator  to  15°  North  and  from  5° 
to  20°  West.  The  area  of  interest  included  the  cities  of  Monrovia  (Liberia), 
Freetown  (Sierra  Leone),  and  the  entire  Liberian  coast.  At  a  resolution  of  1/5  of 
a  degree,  the  model  was  run  on  76  by  76  grid  points.  The  boundary  conditions 
were  setup  using  the  database  of  global  results  available  via  the  GUI.  In  all, 
35  tidal  stations  were  found  in  this  domain,  but  only  22  were  retained.  The  others 
were  excluded  based  on  their  location  (protected  harbors,  etc).  The  bathymetry 
was  also  setup  using  the  GUI  and  the  ET0P05  database.  The  barotropic 
(external)  time  step  was  chosen  as  30  seconds  and  the  model  was  run  with  and 
without  wind  forcing.  Time  series  plots  for  a  station  near  Monrovia  are  shown  in 
Figure  6.  Figure  7  shows  sample  streak  plots  for  the  resultant  barotropic 
currents  with  wind  forcing  for  the  coarse  grid  domain.  The  high  resolution 
nested  grid  extended  from  5°  to  8.4°  North  and  8.7°  to  13.2°  West.  The  grid 
resolution  was  four  times  finer  than  the  coarse  grid  or  1/20  of  a  degree.  The 
selected  domain  yielded  a  89  by  69  grid  with  two  open  boundaries  (West  and 
South).  The  boundary  conditions  were  provided  from  the  saved  output  of  the 
coarse  model  run  at  the  open  boundaries  of  the  nested  domain.  Data  from  only 
two  tidal  gages  was  available  which  was  assimilated  similarly.  The  barotropic 
time  step  was  reduced  to  8  seconds.  A  refined  coastline  and  an  accurate  high 
resolution  bathymetry  was  used.  A  sample  streak  plot  and  time  series  for 
Monrovia  are  shown  in  Figures  8  and  9  respectively.  As  before,  fine  resolution 
results  capture  greater  detail  and  variability  than  the  coarse  resolution  results. 

FUTURE  RESEARCH: 

In  building  on  our  FY95-96  accomplishments,  future  plans  include: 

Task  1:  Validation  of  the  2D  model  with  Currents.  The  fully  integrated  2D 
model  has  already  been  validated  against  tidal  station  data  and  sea  level 
observations,  and  only  the  barotropic  currents  need  to  be  validated.  The 
barotropic  influences  weigh  heavily  in  the  coastal  and  shallow  waters.  The  model 
predictions  will  be  validated  with  any  available  observations  in  a  suitable  area. 

Task  2:  Application.  Validation,  and  Integration  of  3-D  Baroclinic  Model. 
Continue  the  applications  and  validation  of  the  3-D  baroclinic  model  developed  at 
CU  under  ONR  sponsorship.  This  is  a  sigma-coordinate  based,  comprehensive 
physics  3-D  tidal  model  (with  data  assimilation)  used  to  provide  estimates  of  tidal 
current  structure  (both  vertical  and  horizontal)  in  coastal,  marginal,  and 
semi-enclosed  seas.  The  model  is  currently  running  in  the  Yellow  Sea  at  NRL, 
the  Persian  Gulf  at  NAVOCEANO,  and  is  being  configured  for  the  South  China 
Sea  at  CAST. 
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Figure  6.  Time  Series  plots  for  Monrovia  (  Coarse  Grid  ) 
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Figure  7.  Streak-Plot  of  Tidal  Currents  for  Coarse  Grid  Domain 
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Figure  8.  Streak-Plot  of  Tidal  Currents  for  Fine  Grid  Domain 
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TOTAL  ELEVATION  ABOVE  MEAN  SEA  LEVEL  PRESSURE 
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Figure  9.  Time  Series  plots  for  Monrovia  (  Fine  Grid  ) 
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Tidal  currents  are  important  in  many  applications  including  interpreting  ADCP 
observations.  While  high  resolution  is  preferable,  because  of  the  smaller  scales 
involved  (compared  to  tidal  SSH),  most  often  the  data  needed  for  initialization 
and  forcing  of  high  resolution  tidal  models  are  limited.  There  is  also  a  need  to 
nest  high  resolution  tidal  current  models  of  the  near  coast  to  basin-wide  coarse 
resolution  models  to  provide  the  needed  boundary  conditions.  It  is  with  these  in 
mind  that  CU  will  develop  a  3-D  tidal  current  model  based  on  the  same  approach 
as  the  2-D  tidal  SSH  model,  namely  data  assimilation  to  force  realism  and 
relocatable-nestable  capability.  The  model  will  be  applied  to  coastal, 
semi-enclosed,  and  marginal  seas  of  Navy  interest.  MSU  will  build  the  GUI.  The 
RME  will  be  enhanced  to  support  sigma  coordinates  and  boundary  conditions  in 
the  vertical.  If  reliable  transport  estimates  are  available  they  will  be  prescribed. 
If  not,  a  mean  geostrophic  flow  will  be  used.  The  model  will  be  validated  with 
current  meter  observations  where  available,  and  nested  high  resolution  studies 
would  be  conducted  if  suitable  data  are  available. 

Task  3:  Validation  and  Comparison  of  ADCIRC  Model.  NAVOCEANO's 
ADCIRC  Finite  Element  Model  output  will  be  integrated  into  the  RME  for  the 
purposes  of  model  evaluation  and  comparison  against  available  observational  data 
and  other  appropriate  criteria.  Validation  and  inter-model  comparison  between 
the  CU  model  and  ADCIRC  will  be  performed  in  close  collaboration  with 
NAVOCEANO  and  as  such,  will  be  conducted  in  earmarked  regions  of  interest 
and  with  suitable  parameterization  developed  in  collaboration  with 
NAVOCEANO.  General  guidelines  will  also  be  developed  for  the  above,  which 
will  then  be  incorporated  into  the  RME  interface  for  the  future  validation  of 
other  relocatable  models. 


18 


OTHER  RESEARCH  INTEREST;  MODELING  WITH  DATA 
ASSIMILATION  IN  THE  NORTH  ATLANTIC  (DAMEE) 

by  Dr.  Avichal  Mehra 


SCIENTIFIC  ADVISOR;  Dr.  David  E.  Dietrich 
OBJECTIVE; 

To  show  the  strengths  and  weaknesses  of  different  modeling/assimilation 
approaches. 

APPROACH; 

For  the  past  two  years,  DAMEE-NAB  has  addressed  the  climatological  behavior 
of  different  prognostic  circulation  models  as  per  a  fixed  list  of  well  documented 
properties  of  the  North  Atlantic  Basin.  Some  forecast  experiments  on  the  target 
domain  of  the  North  Atlantic  Subtropical  Gyre  have  also  been  performed.  These 
experiments  entailed  various  data  assimilation  methods  and  emphasized  synoptic 
time-scale  predictions.  All  such  experiments  have  been  performed  at  a  low 
resolution  of  1/2°  with  fewer  than  20  vertical  levels.  Future  participants  of 
DAMEE-NAB  will  address  the  sensitivity  of  mesoscale  forecasts  to  variations  in 
climatological  measures,  and  also  establish  basin-scale  predictive  capabilities  of 
general  circulation  models  relative  to  persistence  and  climatology  at  higher 
resolutions  (~1/10°)  horizontally  and  a  larger  number  of  vertical  levels  (up  to 
40).  Some  other  goals  include  identifying  the  strengths  and  weaknesses  of 
different  classes  of  numerical  models  and  various  data  assimilation  techniques. 
Exploring  and  implementing  the  best  data  assimilation  methods  with  coupled 
models  will  be  the  final  step  for  accurate,  efficient  forecasts  of  the  North  Atlantic 
Basin. 

Recently  DieCAST,  with  a  modified  Arakawa  "a"  grid  scheme,  was  successfully 
implemented  in  the  North  Atlantic  Basin  at  one-half  degree  resolution  with  ten 
vertical  levels.  The  domain  extended  from  15°S  to  65°N,  and  from  the  equator  to 
80°W.  It  included  the  entire  North  Atlantic  basin  and  Atlantic  Equatorial  Region, 
but  excluded  the  Gulf  of  Mexico.  The  vertical  viscosity  was  determined  by  a 
method  based  on  Richardson  number,  while  realistically  low  horizontal 
diffusivities  were  specified  to  model  the  circulation.  Results  indicate  active 
inertial  eddy  fields  caused  by  implementation  of  higher-order  interpolations  and 
modifying  the  incompressibility  algorithm  to  substantially  reduce  dispersion. 
Even  at  a  coarse  resolution  of  one-half  degree,  the  Gulf-Stream  separation  seems 
robust  (see  Figure  10)  and  other  well-known  features  of  the  North  Atlantic 
Subtropical  Gyre  were  reasonably  simulated.  Particularly  striking  are  realistic 
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Equivalent  Hydrostatic  SSH  (cm)  at  Day  6820 
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Fiaure  10 


detailed  features  associated  with  the  Flemish  Cap  (Figure  10).  This  includes  an 
anticyclonic  Taylor  column;  an  anticyclonic  Gulf  Stream  water  penetration  east 
of  the  Flemish  cap  into  the  Labrador  Sea,  between  two  semipermanent  small-scale 
cyclonic  features;  and  a  southward  moving  central  Labrador  Current  branch  in 
the  Flemish  Channel  on  the  west  side  of  the  Flemish  Cap,  with  the  full  observed 
three  branch  system  clearly  shown  in  the  Labrador  Sea.  This  research  involves 
the  following  three  tasks. 

1.  Comparison  with  Observations.  DieCAST  will  be  implemented  at 
1/10°  resolution  with  40  z-levels  in  the  North  Atlantic  Basin  including  the  Gulf 
of  Mexico.  As  noted  above,  reasonable  results  have  already  been  obtained 
using  1/2°  resolution  and  10  z-levels.  To  be  sure  that  the  model  is  converged, 
an  intermediate  resolution  run  will  also  be  performed.  The  NAVOCEANO 
GDEM  climatology  will  be  used  along  with  ECMWF  (1993)  winds  with 
variability.  Some  results  that  will  be  analyzed  include  the  mean  and  seasonal 
cycle  of  Florida  Straits  water  at  27.5°N,  eddy  kinetic  energy  distributions  at 
various  depths,  mean  and  variance  of  Gulf-stream  path,  and  separation  at  Cape 
Hatteras.  Certain  model  vertical  cross-sections  at  27°N  and  55°W  will  also  be 
compared  and  contrasted  against  observations.  Model  outputs  adequate  to 
show  performance  criteria  will  be  included.  In  addition,  because  transient 
eddies  may  be  critical  to  the  maintenance  of  the  Labrador  Current  between  the 
Grand  Ba^s  and  Cape  Hatteras  and  thus  also  to  the  realistic  GS  separation, 
diagnostics  to  assess  the  roles  of  transient  eddies  ("negative  viscosity"  alias 
"Neptune  effect")  will  be  included. 

2.  Study  Thermohaline  Circulation.  The  Labrador  Sea  is  one  of  the  most 
intensive  water  mass  transformation  regions  of  the  world.  It  affects  the 
thermohaline  circulation  of  the  entire  NAB.  In  the  DAMEE  experiment,  this 
is  treated  through  specification  of  a  northern  buffer  zone  climatology  to  obtain 
proper  deep  water  inflows  from  the  Labrador  Sea.  To  model  the  water  mass 
transformations  more  directly,  the  domain  may  have  to  be  extended  to  75°N 
(if  resources  allow).  The  resulting  thermohaline  effects  will  be  addressed  by 
(a)  exploring  the  sensitivity  of  the  circulation  to  turbulent  vertical  diffusion, 
downward  diffusion  of  latitudinally  differential  density  from  surface  heat  and 
salinity  fluxes  creates  available  potential  energy  for  the  thermohaline 
circulation;  (b)  maintaining  a  diagnostic  time  mean  of  model  surface  heat  and 
salinity  fluxes  at  all  horizontal  grid  points;  (c)  adding  the  long  term  mean 
model  heat  and  salinity  heat  fluxes  to  the  surface  layer  each  time  step;  (d)  after 
adding  the  long  term  mean  fluxes,  restoring  to  GDEM  with  O  (180)  day 
restoring  time  scale  each  step;  and  (e)  using  the  streamfunction  for  the  zonally 
averaged  flow  to  diagnose  thermohaline  effects  on  the  meridional  circulation 

3.  Implement  Data  Assimilation.  The  DieCAST  ocean  model  capabilities  will  be 
enhanced  by  incorporating  the  best/most  practicable  data  assimilation. 
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OTHER  RESEARCH  INTEREST:  A  NESTED  HYBRID 
MODELING  SYSTEM 
by  Dr.  Avichal  Mehra 


OBJECTIVE: 

To  deliver  an  accurate  and  efficient  hybrid  modeling  system  for  nowcast/forecast 
applications  for  ocean/tidal  currents  and  sea  surface  height. 

APPROACH: 

DieCAST  is  an  efficient,  robust,  accurate  model  in  application  to  deep  water  and 
to  coupled  coastal  flows,  as  well  as  to  island  wakes  and  detailed  high  resolution 
flows  in  a  fjord.  These  positive  attributes  result  from  its  rigid-lid  approximation, 
z-level  vertical  coordinate,  conservative  control  volume  numerics,  and  special 
accurate  treatment  of  Coriolis  and  internal  wave  propagation  terms.  A  free 
surface  with  explicit  time-split  barotropic  mode  has  been  recently  added. 

CURReNTSS  is  a  sigma  coordinate,  finite  difference,  explicit,  vertically- 
integrated  barotropic  model  assimilating  tidal  data  from  coastal  gages  and  from 
bottom  pressure  gages.  It  is  fully  non-linear  and  incorporates  direct  astronomical 
forcing,  and  can  also  include  surface  forcing  from  atmospheric  models  to  predict 
storm  surges.  The  model  can  be  run  with  any  number  of  components  including 
long  term  and  compound  tides.  It  employs  a  simple  data  assimilation  procedure 
by  replacing  the  model  predicted  SSH  at  pre-determined  intervals  by  a  weighted 
sum  of  the  model  prediction  and  the  observed  SSH  from  the  tide/bottom  pressure 
at  that  point,  the  weights  determined  a  priori.  The  model  is  well  formulated, 
carefully  calibrated  and  adequately  validated  for  a  number  of  regions.  The  basic 
problem  is  that  while  bottom  boundary  fit  approaches  (such  as  sigma  coordinates) 
may  be  attractive  for  very  shallow  regions,  they  are  not  designed  for  accurate 
representation  of  baroclinic  pressure  gradient  and  vortex  stretching  terms,  which 
involve  gradients  of  density  and  horizontal  velocity  in  regions  of  large  amplitude 
and  steep  topography. 

Dr.  Paul  Martin  (NRL)  has  noted  the  stair-step  approximation  can  cause 
spurious  local  horizontal  divergence  that  does  not  occur  in  sigma  coordinate 
models.  He  also  noted  that  when  the  pycnocline  appears  as  a  sharp  "horizontal" 
front  because  of  sloping  sigma  surfaces,  horizontal  advection  errors  can  be 
amplified.  This  can  lead  to  large  spurious  oscillations  when  using  realistic 
diffusivities.  However,  large  diffusivities  are  undesirable  because  they  can  cause 
unrealistic  cross-pycnocline  mixing.  In  z-level  models,  horizontal  advection 
errors  are  smaller,  because  the  pycnocline  is  quasi-horizontal  and  thus  easier  to 
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resolve  horizontally,  while  vertical  advection  errors  are  small  simply  because 
buoyancy  suppresses  vertical  velocity  properly  near  thermoclines  in  z-level 
models.  It  thus  seems  appropriate  to  blend  sigma  coordinate  coastal  models  with 
deep  water  models.  This  is  presently  being  done  by  Dr.  D.S.  Ko  in  coupling  the 
Princeton  Model  to  the  Navy  Layered  Model  and  by  Dr.  Martin  within  a  single 
model.  DieCAST  provides  another  option,  because  it  can  be  used  with  the  high 
vertical  resolution  needed  for  accurate  coastal  modeling  and  can  conveniently 
address  air-sea  exchanges.  The  goal  is  to  use  DieCAST  to  provide  deep  water 
open  boundary  conditions  to  the  coastal  model  and/or  use  the  models  as  a  fully 
coupled  system.  Tasks  involved  are: 

1.  Blend  of  sigma  and  z-level  co-ordinates.  Techniques  to  modify  sigma 
coordinate  surfaces  instead  of  the  often  used  topographic  filter  approach  will 
be  investigated  to  improve  pressure  gradient  and  vortex  stretching  numerics. 
Such  modifications  reduce  vertical  coordinate  surface  slopes  in  deep  water, 
without  affecting  bottom  topography.  This  may  avoid  the  need  for  topographic 
filters.  Two  approaches  will  be  considered.  The  first  approach  (filters)  will 
start  with  conventional  sigma  coordinate  surfaces.  These  will  be  filtered  by 
iterative  application  of  a  Laplacian  filter  designed  to  decrease  large  coordinate 
surface  slopes,  but  only  when  they  are  away  from  the  bottom;  the  slope 
damping  will  be  proportional  to  the  distance  from  the  bottom.  The  second 
approach  is  to  define  each  individual  vertical  coordinate  surface  as  a  function 
of  depth  such  that  it  is  quasi-horizontal  in  areas  deeper  than  a  specified  depth. 
Thus,  its  depth  function  will  asymptotically  approach  a  specified  value  as  the 
depth  increases.  The  resulting  vertical  surfaces  will  substantially  deviate  from 
sigma  coordinates,  but  little  or  no  change  will  be  needed  to  use  their  metrics  in 
conventional  sigma  coordinate  models. 

2.  Flux  conserving  nesting.  The  first  requirement  by  a  nested  submodel  is  the 
normal  boundary  velocity,  which  must  have  integrated  zero  volume  flux.  The 
other  fundamental  boundary  quantities  that  must  be  conserved  are  momentum, 
salinity  and  heat  fluxes.  It  may  be  simplest  and  best  to  specify  the  actual  nested 
fluxes  of  momentum,  salinity  and  heat  from  the  z-level  coarse  model.  This 
eliminates  guess  work  for  open  boundary  flux  specification  and  allows  a 
consistent  interior  nested  domain  solution  to  develop.  Initially,  the  coupling 
will  be  implemented  between  the  Arakawa  "c"  grid  version  of  DieCAST  and 
CURReNTS.  Later,  the  Arakawa  "a"  grid  version  will  be  used.  Such  a  hybrid 
system  would  take  advantage  of  the  strong  features  of  both  z  and  sigma  level 
modeling  techniques.  It  would  provide  less  dispersive,  more  accurate  and 
efficient  results  for  deep  water,  and  their  use  for  accurate  predictions  for 
coastal  currents,  tides  and  sea  surface  elevations  from  sigma-level  model. 
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A  STUDY  OF  EDDY  DYNAMICS  IN  THE  NAVY  LAYERED 
PACIFIC  OCEAN  MODEL  DATA 
by  Dr.  Zhifan  Zhu 

SCIENTIFIC  ADVISOR:  Dr.  Robert  Moorhead,  MSU-NSF  Engineering 

Research  Center 

OBJECTIVE: 

Eddies  are  prominent  ocean  features  and  play  a  critical  role  in  ocean 
circulation  modeling  and  prediction.They  are  time-dependent  oceanic  phenomena 
that  exist  ubiquitously  in  ocean  circulations  and  are  often  associated  with  strong 
barometric  and  thermal  variabilities  such  as  low/high  pressure  centers  and  radial 
temperature  gradients.  Various  interactions  and  coupling  effects  between 
mesoscale  eddies,  ocean  currents,  and  heat  transport  have  long  been  studied. 
However,  due  to  the  complex  mechanisms  involved  in  eddy  dynamics,  many  eddy 
aspects  are  not  understood;  how  eddies  contribute  to  and  interact  with  the 
dynamics  of  other  oceanic  features  remains  largely  unresolved.  Characterizing 
eddy  dynamics  over  time  and  space  has  thus  been  of  keen  interest  to  physical 
oceanographers  in  recent  decades  and  continues  to  be  a  major  research  topic  in 
physical  oceanography 

Modeling  of  global  scale  ocean  circulations  has  recently  become  an  important 
means  of  oceanographic  studies.  In  particular,  models  with  eddy-resolving 
capability,  have  played  a  central  role  in  global  ocean  prediction.  The  time- 
varying  three-dimensional  model  outputs  consist  of  high  resolution  data  samples 
of  the  ocean  circulation  fields  and  are  well  suited  to  the  study  of  eddy  dynamics. 
Here  the  objective  was  to  conduct  an  investigation  of  eddies  in  two  regions  of  the 
Navy  Layered  Pacific  Model  data  from  1983  to  1987. 

APPROACH: 

The  approach  is  based  on  an  eddy  extraction  algorithm  developed  for 
extracting  and  visualizing  eddies  from  time-dependent  flow  fields.  The  Navy 
Layered  Ocean  Model  was  developed  by  the  Naval  Research  Laboratories  (NRL) 
at  Stennis  Space  Center.  The  data  used  was  fi-om  a  1/8°  eddy-resolving  model  of 
the  Pacific  Ocean  north  of  20°S.  Artificial  solid  boundaries  are  placed  at  20°S 
and  in  the  straits  at  the  southern  boundary  of  the  Indonesian  Archipelago. 
Horizontally,  the  model  covers  the  area  from  110°E  to  78°W  and  from  20°S  to 
62°N,  with  a  horizontal  resolution  of  0.125°  x  0.176°  (lat.  Ion)  for  each  variable. 
Vertically,  it  has  6  layers  with  mean  interface  depths  of  135,  320,  550,  800,  and 
1050  meters.  The  other  major  features  of  the  model  include  nonlinear  primitive 
equations,  a  free  surface,  a  semi-implicit  time  scheme,  full-scale  bottom 
topography  in  the  lowest  layer,  and  an  arbitrary  coastline  geometry. 
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Starting  from  rest,  the  model  was  spun  up  to  statistical  equilibrium  at  1/4° 
resolution  using  the  Hellerman-Rosenstein  monthly  wind-stress  climatology,  and 
then  continued  another  15  years  at  1/8°  resolution.  Finally,  it  was  run  for  1981- 
1989  using  monthly  averaged  ECMWF  1000-mb  winds  with  the  1981-1989  mean 
replaced  by  the  annual  mean  from  Hellerman-Rosenstein.  No  oceanic  data  were 
assimilated  in  this  purely  wind-forced  simulation.  Verification  of  the  model  was 
conducted  using  satellite  altimeter  data  from  the  Geosat  Exact  Repeat  Mission 
(ERM).  Dr.  James  Mitchell  (NRL)  compared  the  statistics  of  the  simulated  and 
observed  SSH  variations  for  the  Kuroshio  Extension;  while  Gregg  Jacobs  (NRL) 
analyzed  the  two  SSH  variabilities  through  frequency-wavenumber 
decomposition.  Both  studies  reported  good  correlation  between  model  and 
altimeter  SSH  variations. 

The  ou^ut  consists  of  layer  thickness  deviation  '/f  and  horizontal  transport 
velocity  {V  =  (U,V)}  for  a  3.05  model  day  interval  at  each  layer  and  horizontal 
location.  The  relationship  of  flow  velocity  and  transport  velocity  Vj,  at  the  kth 
layer  is  given  by: 

Vk  =  hkVk 


where  hk  is  the  layer  thickness  which  can  be  derived  from  Hk  and  the  reference 
depth.  The  total  data  size  for  a  four-year  run  is  about  22  gigabytes. 

Two  regions  were  selected  for  the  study.  Region  A  is  from  110°W  to  145°W 
and  from  15°N  to  62°N  along  the  west  coast  of  North  America,  and  Region  B 
extends  from  120°E  to  180°E  and  from  20°N  to  55°N,  containing  the  Kuroshio 
Current  near  Japan. 

The  algorithm  extracts  eddies  from  the  time-varying  volumetric  flow  fields 
and  stores  them  in  an  abstract  form.  Each  extracted  eddy  consists  of  a  sequence 
of  eddy  instances  (in  a  3D  spatial  sense)  from  emergence  to  decay.  Each  instance 
is  further  decomposed  into  a  stack  of  horizontal  slices  (2D  layers).  The  number 
of  slices  can  vary  from  one  to  six.  A  minimum  life-span  threshold  was  imposed 
in  this  calculation;  eddies  that  existed  less  than  two  months  were  excluded. 

RESULTS: 

Given  a  two-month  life  time  criterion,  207  eddies  in  Region  A  and  314  eddies 
m  Region  B  were  extracted  from  the  four  year  data.  The  most  persistent  eddy  in 
Region  A  had  a  life-span  of  311.1  days  while  in  Region  B  it  was  741.2  days. 
These  521  eddies  form  the  basis  of  the  following  analysis. 
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a.  Trajectories.  Figures  la  and  lb  show  the  eddy  trajectories  visualized  by  tracing 
the  centers  of  the  tracked  eddies.  The  red  represents  cyclonic  eddies  and  the  yellow 
anticyclonic  eddies.  The  grayish  background  is  the  bathymetry.  Notice  that  the 
images  show  the  superimposing  eddy  migration  paths  during  the  four  year  period. 
The  time  information  is  lacking.  This  information  can  be  effectively  investigated 
by  an  interactive  animation  process,  which  is  not  permitted  here. 

The  eddies  in  Region  A  are  concentrated  off  the  California  Coast  with  the  densest 
areas  below  40°N.  While  below  30°N  the  eddy's  movements  are  consistently  in  the 
horizontal  direction,  the  eddies  from  30°N  to  40°N  near  the  coast  moved  in 
multiple  directions.  Long  migration  paths  are  easily  seen  away  from  the  coast 
below  30“N.  These  eddies  were  created  from  the  water  flowing  parallel  to  the 
North  Equatorial  Current.  The  existence  of  alternating  cyclonic  and  anticyclonic 
eddies  from  15*N  to  35“N  (west  of  130’W)  suggests  significant  horizontal  currents 
in  that  area,  causing  the  eddies  to  be  pinched  off  above  and  below  the  currents. 
The  eddies  at  the  high  latitudes  (above  40°N)  were  mostly  anticyclonic,  which  tend 
to  meander  in  the  littoral  waters.  No  (long-term)  eddies  were  found  above  35'N. 

The  average  migration  distance  of  the  eddies  in  Region  A  is  454.46  kms  and  the 
mean  speed  is  5.07  cms/s.  In  general,  the  eddies  at  the  lower  latitudes  showed 
larger  speeds.  The  speed  difference  between  the  cyclonic  and  the  anticyclonic 
eddies  is  trivial,  although  there  were  twice  as  many  cyclonic  as  anticyclonic  eddies. 

Region  B  can  be  divided  into  three  separate  areas:  Sea  of  Japan,  Sea  of  Okhotsk, 
and  the  NW  Pacific  region.  Since  higher  resolution  models  of  the  Sea  of  Japan 
exist,  and  the  Sea  of  Okhotsk  is  only  partially  covered,  the  analysis  here  is  focused 
on  the  Pacific  region  which  includes  the  Kuroshio  Extension.  The  eddy  trajectories 
are  in  an  elongated  area  around  the  Kuroshio  Current.  The  number  of  cyclonic 
and  anticyclonic  eddies  is  almost  equal.  The  highest  eddy  density  occurred  in  the 
area  to  the  west  of  the  Shatsky  Rise  and  extended  northeastward  to  the  Emperor 
Seamount  Chain  (at  ~170‘’E).  The  trajectories  are  non-linear,  revealing  the  great 
energy  and  complexity  of  currents  surrounding  the  Kursohio  Extension.  Beyond 
the  Emperor  Seamount  Chain,  most  eddies  wandered  in  small  areas  creating 
twisted  trajectories.  Most  of  the  coastal  eddies  can  be  found  off  the  Tsugaru  Strait. 
The  tangled  lines  show  the  superimposing  eddy  trajectories  over  the  four  years. 
The  long  and  stretched  pathlines  are  a  distinct  pattern  in  the  western  part  of  the 
Kuroshio  (below  ~3rN),  particularly  across  the  Izu  Ridge.  These  eddies  were 
created  from  rings  pinched  off  the  horizontal  stream.  The  average  distance  of 
eddy  migration  in  this  region  is  434.50  kms,  and  the  average  speed  is  4.80  cms/s. 
Noteworthy  is  that  Aoki  reported  a  typical  speed  of  4  cm/s  for  two  strong  cyclonic 
eddies  tracked  in  the  first  year  (from  November  1986  to  November  1987)  Geosat 
ERM  altimeter  data  around  the  Kuroshio  Extension. 
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Figure  la.  Region  A. 
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b.  Time  sequences  of  eddv  densities.  It  is  known  that  eddies  carry  significant 
amounts  of  energy  in  various  interactions  with  the  oceanic  and  atmospheric 
fields.  The  eddy  field  has  a  direct  influence  on  the  global  heat  budget  and  surface 
and  subsurface  variabilities,  and  vice  versa.  Such  a  relationship  can  be  examined 
by  observing  the  density  of  the  eddy  fields  over  time.  The  density  of  the  eddy 
field  is  represented  by  the  number  of  eddies  that  exist  at  a  particular  time.  In 
Figure  2  are  two  curves  showing  the  results  for  each  model  day  during  the 
four  years.  The  eddy  field  in  Region  B  is  generally  correlated  with  the  season, 
e.g.,  more  eddies  exist  in  the  summer  time  than  during  the  winter.  The  annual 
peaks  are  around  August-September.  A  small  peak  in  January  1985  indicates  an 
'anomaly'.  The  variations  are  generally  in  agreement  with  the  seasonal  cycle  as 
estimated  by  Tai  and  White  (1990)  using  the  rms  sea  level  variability  from 
GEOSAT.  Since  the  eddy  field  around  the  Kuroshio  Extension  is  primarily 
forced  by  the  instabilities  of  the  strong  current,  one  would  expect  a  similar  cycle 
in  the  current.  The  eddy  field  in  Region  A  shows  a  seasonal  correlation  only  in 
the  last  two  years  (1985-86).  In  particular,  a  sudden  drop  in  eddy  density 
occurred  in  the  summer  of  1983  and  caused  an  180°  out  of  phase  until  the  spring 
of  1985.  This  phenomenon  may  be  related  to  the  strong  El  Nino  event  ending  in 
1983,  the  fact  that  low  latitude  variability  is  less  seasonal  by  itself  (Pond  and 
Richard,  1983),  and  the  fact  that  the  amount  of  vertical  motion  in  the  littoral 
areas  may  mask  the  seasonality. 

The  auto-correlation  functions  of  the  two  time  sequences  in  Figure  2  are  shown 
in  Figure  3.  The  alternating  pattern  for  Region  B  is  a  strong  indication  of  the 
eddy  field  varying  in  seasonal  cycles. 

c.  Sea  surface  height  (SSHJ  of  eddv  centers.  The  SSHs  are  computed  at  the  top- 
layer  eddy  centers,  and  are  averaged  separately  over  the  cyclonic  and  anticyclonic 
eddies  at  each  time  step.  They  indicate  the  variability  of  the  sea  surface  elevation 
of  the  eddy  field.  Figure  4  shows  the  variations  of  SSHs  over  time.  In  Region  A, 
the  SSHs  of  the  cyclonic  eddies  and  anticyclonic  eddies  vary  in  similar  patterns 
over  time,  indicating  that  the  SSHs  of  the  eddy  field  closely  follow  the  overall 
SSHs  in  the  domain.  This  is  probably  because  the  amplitude  of  sea  surface 
variations  at  the  eddy  centers  caused  by  the  vertical  water  transportation  is  only  a 
fraction  of  the  domain  SSH  deviations.  In  Region  B,  the  correlations  between  the 
cyclonic  and  anticyclonic  SSH  patterns  are  not  as  obvious  as  in  Region  A,  and  the 
magnitude  of  the  SSH  variations  is  about  twice  as  large  as  that  in  Region  A,  thus 
indicating  the  model  output  being  consistent  with  climatology  because  eddies  are 
more  energetic  towards  Ae  west  in  the  North  Pacific  Ocean. 
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eddies 


Figure  2.  Time  sequences  of  the 
number  of  tracked  eddies. 


Figure  3.  Auto-correlation  function 
of  the  time  sequences  of  the  number 
of  tracked  eddies. 


Figure  4.  Sea-surface  height  at  the  eddy  centers. 


FUTURE  RESEARCH: 

Future  efforts  will  be  directed  at  creating  a  system  that  can  intelligently 
digest  large  amount  of  raw  data  and  save  the  most  relevant  information  as 
features  in  an  object-oriented  database  for  various  oceanographic  applications. 
The  system  should  be  able  to  provide  a  high-level  description  of  ocean  physics  at 
a  given  time  and  place.  Such  information,  based  on  the  features  stored  in  the 
database,  should  be  sufficient  for  a  user  to  characterize  the  physical  domain  under 
study. 
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ENVIRONMENTAL  DATA  VISUALIZATION  SYSTEM 
VERSION  2.0  (ENVIS  2.0) 
by  Dr.  Zhifan  Zhu 

SCIENTIFIC  ADVISOR:  Dr.  Harsh  Anand,  Center  for  Air  Sea  Technology 
TECHNICAL  ADVISOR:  Mr.  James  Corbin,  Center  for  Air  Sea  Technology 

OBJECTIVE: 

The  development  of  ENVIS  2.0  is  based  on  ENVIS  1.0  and  the  Unified  Air 
Sea  Visualization  System  (UASVS).  ENVIS  was  developed  as  a  tool  that 
emphasized  feature  detection  capabilities,  and  it  received  positive  feedback  from 
oceanographic  data  visualization  applications.  However,  the  software  was  built 
upon  Silicon  Graphics  workstations  and  its  data  interface  capability  was  limited. 
UASVS  focused  on  the  gridding  structures  present  in  different  ocean  model  data 
sets.  It  allows  one  to  map  the  original  grids  to  a  single,  unified  grid  so  users  can 
easily  visualize  and  compare  the  results  of  various  models.  Although  both 
ENVIS  and  UASVS  possess  attractive  functionalities,  their  non-objected-oriented 
nature  makes  them  difficult  to  upgrade. 

The  objective  of  ENVIS  2.0  is  to  incorporate,  refine  and  expand  the  data 
visualization  capabilities  of  ENVIS  1.0  and  UASVS,  using  object-oriented 
designs.  That  will  provide  flexibility  in  maintaining  and  updating  the  system. 
The  immediate  improvements  include  simultaneously  visualizing  time-varying  3D 
data  objects  with  different  time  and  spatial  scales,  without  interpolating  data,  i.e., 
data  objects  will  be  kept  in  their  original  resolutions  and  formats;  interactive 
edition  of  graphic  object  attributes,  e.g.,  individual  colormap  and  texture  will  be 
associated  with  each  graphic  object;  interactive  edition  of  text  labels,  i.e., 
allowing  the  user  to  generate  and  edit  text  labels  in  the  fly;  animation  of  graphic 
objects  in  time  and  space  and  the  time-space  animation  information  can  be  saved 
and  loaded  for  successive  rendering;  and  allowing  output  to  rendered  images  in 
various  formats  to  generate  postscript  visual  presentation.  The  ultimate  goal  of 
ENVIS  is  to  make  it  a  scalable  visualization  structure  such  that  various 
applications  can  be  built  effectively  and  efficiently.  Applications  may  include 
stand-alone  visualization  tools  on  different  platforms,  and  functional  modules 
integrated  in  a  distributed  modeling  and  data  management  environment  using 
CORBA/EDGE  technology  currently  being  developed  at  CAST. 

APPROACH: 

In  order  to  achieve  interoperability  and  portability  using  object-oriented 
techniques,  the  essential  functionalities  of  ENVIS  2.0  are  grouped  into  hierarchies 
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of  classes.  These  hierarchies  include  data  object,  graphic  object,  render  class, 
GUI,  animation,  image,  and  label  classes.  Since  we  are  dealing  with 
heterogeneous  data  sources  (e.g.,  grids  and  formats),  a  data  class  hierarchy  needs 
to  be  developed.  The  root  class  is  an  abstract  class  under  which  specific  data 
interface  classes  are  created.  These  data  classes  have  unified  API  with  other  class 
groups  and  are  able  to  handle  various  data  sources  implicitly, 

The  graphic  object  classes  consist  of  polygons,  mapping  routines  to  map  the 
data  objects  into  graphic  objects,  and  rendering  routines.  These  routines  are 
primarily  written  in  OpenGL,  a  predominant  graphic  language  supported  by  most 
graphic  workstations.  Render  classes  consist  of  a  pool  of  machine  dependent 
renders.  First,  the  rendering  classes  are  designed  for  each  type  of  machine,  e.g., 
SGIs  and  SUNs.  Then  sub-classes  are  developed  with  each  parent  class.  For 
SGIs,  the  sub-classes  will  be  based  on  their  graphic  hardware  capabilities,  e.g.. 
Indy  8-bit,  Indy  24-bit,  Indigo2  Extreme,  Indigo2  High  Impact,  etc.,  to  assure  an 
economical  fit.  Under  the  same  hardware  render  class,  2D  and  3D  render  classes 
are  created  for  different  applications. 

The  development  of  GUI  classes  is  not  a  major  task  in  ENVIS  2.0.  The 
effort  is  leveraged  with  other  projects  at  CAST  to  create  a  unified  GUI  library. 
In  particular,  ENVIS  2.0  can  be  integrated  as  a  visualization  model  in  an 
application,  such  as  Weather  Watch.  No  top  level  GUI  is  needed.  Animation, 
image,  and  label  classes  are  supporting  classes  to  facilitate  the  rendering  process 
and  creation  of  technical  presentations. 

RESULTS: 

A  prototype  of  ENVIS  2.0  has  been  developed.  It  consists  of  libraries  of 
classes  grouped  in  functional  hierarchies.  Additions  of  new  classes  within  a  class 
hierarchy  do  not  affect  the  interface  with  other  class  hierarchies.  Figure  5 
displays  the  relationship  among  these  class  hierarchies. 


Figure  5.  Class  relationships. 
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Corresponding  GUIs  have  been  developed  to  integrate  these  functional  objects 
into  a  visualization  tool  that  the  user  can  steer.  Figure  6  shows  ENVIS  2.0  GUIs 
which  are  implemented  with  a  Motif  class  library.  The  mainwindow  contains  a 
SGI  render  and  toplevel  menu  bar  to  launch  other  GUIs.  DataManager  lets  the 
user  manage  data  objects,  e.g.,  loading  and  deleting.  It  also  allows  selection  of 
section  of  data  from  a  data  object  for  mapping.  ViewManager  controls  the 
mapped  graphic  objects,  e.g.,  selecting  or  deselecting  graphic  objects,  or 
modifying  the  graphic  attributes  such  as  colormap  and  opacity  of  selected  graphic 
objects.  Graphic Attr  maintains  a  workspace  to  allow  the  user  to  modify  graphic 
attributes.  LabelManager  controls  the  creation,  deletion  and  modification  of  label 
objects  such  as  string  labels,  time  information  display,  and  colorbars. 
ImageManager  allows  one  to  save  currently  rendered  images  in  any  supported 
image  formats.  AnimateManager  provides  controls  of  the  animation  in  both  time 
and  space.  For  temporal  animation,  the  user  can  select  any  time  window  and  time 
increment  to  animate  the  timed  graphic  objects  in  both  forward  and  backward 
directions.  For  space  animation,  the  user  can  pre-register  viewing  angles  and 
positions  interactively  and  play  back.  The  time  and  space  animation  can  be  used 
simultaneously.  Animation  information,  such  as  time  stamps  and  viewing  angles, 
can  be  saved  to  and  loaded  from  a  simple  text  file  for  session  recording. 


Figure  6.  Snapshot  of  ENVIS  2.0  GUIs. 
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The  Figure  7  shows  an  example  of  using  Envis  2.0  to  visualize  Persian  Gulf 
temperature  and  current  fields  from  the  Princeton  Ocean  Model.  This  version  of 
the  model  has  a  horizontal  curvilinear  and  vertical  sigma  gridding  structure.  The 
temperature  is  visualized  by  color-fill  contours  with  the  color  map  displayed  on 
the  right  side  of  the  image.  The  current  is  rendered  by  streamlines  superimposed 
on  the  temperature  field.  The  terrain  and  bathymetry  are  shown  in  a  different 
color  map.  All  the  rendered  graphics  are  treated  as  graphic  objects.  They  can  be 
'turned  on  and  off  by  selecting  and  deselecting  their  icons  in  the  GUI.  Their 
graphic  attributes  can  be  interactively  modified.  For  example,  the  temperature 
color  contours  and  the  terrain-bathymetry  have  their  own  color  maps  to  enhance 
the  visual  effectiveness.  The  3D  graphic  objects  can  be  easily  transformed 
(translated,  rotated,  and  zoomed  in  and  out)  with  the  mouse.  The  labels 
(including  the  color  map)  can  be  picked  and  placed  at  any  location  in  the  window. 
Their  fonts  can  be  altered  interactively. 

Because  of  the  object-oriented  data  class  design,  different  types  of  data  grids 
and  formats  can  be  easily  supported  and  visualized  in  the  same  volume  of  time- 
space.  There  is  no  centralized  gridding  system.  Data  objects  are  mapped  into 
graphic  objects  directly  from  their  original  grids  and  resolutions.  This  greatly 
reduces  the  cost  of  copying  and  re-formating  data. 

FUTURE  RESEARCH: 

Further  efforts  in  developing  ENVIS  will  be  focused  on  developing  and 
integrating  function  classes  to  be  portable  across  many  display  environments.  To 
do  so,  a  common  interface  approach  will  be  adopted.  This  effort  will  be 
performed  in  conjunction  with  the  Ext;  nsible  Dynamic  GUI  Environment 
(EDGE),  a  research  effort  currently  conducted  by  CAST.  With  this  technique, 
ENVIS  2.0  visualization  capability  will  be  easily  integrated  into  any  modeling  and 
information  management  system  operating  on  most,  if  not  all,  platforms. 
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Figure  7.  Visualization  of  Persian  Gulf  POM  temperature  and  current  using 


OTHER  RESEARCH  INTEREST:  DISTRIBUTED  MULTI¬ 
RESOLUTION  VISUALIZATION 
by  Dr.  Zhifan  Zhu 


OBJECTIVE: 

Creating  better  models  and  obtaining  more  and  better  data  has  always  been  a 
challenge  in  ocean  modeling  and  simulation.  Perceiving,  understanding,  and 
communicating  the  knowledge  contained  in  the  model's  output  or  the  measured 
data  is  critical.  The  fundamental  problem,  as  models  get  larger  and  computers 
get  faster,  is  the  need  to  examine  and  understand  more  data,  in  increasing  detail 
but  in  the  same  amount  of  time.  Distributed  visualization  provides  collaboration, 
concurrentness,  and  computer/visualization  resource  efficiency.  It  is  an  important 
emerging  technology  that  will  enable  scientific  programmers  to  extract 
meaningful  information  rapidly  from  detailed  simulations,  by  eliminating  much 
of  the  disk-to-disk  copying  that  is  now  performed. 

A  typical  situation  of  a  model-visualization  configuration  is  that  the  model  is 
running  on  a  large  and  fast  supercomputer  or  running  in  parallel  over  a  cluster 
of  networked  computers.  Parallel  computation  allows  partition  of  a  large  model 
into  a  set  of  solvers  each  running  on  a  single  processor.  The  results  will  be 
collected  in  real-time  and  sent  to  a  graphic  workstation  for  visualization.  The 
visualization  is  in  fact  a  part  of  model  steering  that  allows  the  modeler  to 
interactively  monitor  the  simulation  by  visualizing  the  results  and  controlling  the 
simulation  according  to  the  visualized  results.  In  this  way,  large  amounts  of  data 
often  need  to  be  transferred  from  the  solvers  to  the  graphic  workstation  through 
the  local  area  network  and/or  wide  area  network.  To  achieve  real  or  near  real¬ 
time  operation,  data  traffic  has  to  be  minimized.  This  can  be  done  with  multi¬ 
resolution  and  data  compression  techniques.  Multi-resolution  means  that  the 
modeler  can  choose  different  resolutions  in  visualizing  the  model  data.  In 
particular,  data  will  be  transmitted  in  a  way  that  visualization  can  begin  at  a 
lower  resolution  and  be  refined  to  higher  resolution.  The  transmission  can  be 
stopped  once  sufficient  information  is  visualized.  It  has  been  identified  that  the 
whole  process  involves  multi-resolution  data  transform  and  compression.  The 
objective  of  this  research  is  to  develop  multi-resolution  and  data  compression 
techniques  that  best  suit  the  distributed  visualization  paradigm. 

APPROACH: 

The  data  transform  plays  two  roles:  packing  the  data  in  a  multi-resolution 
form  that  can  be  progressively  transferred,  and  decorrelating  the  original  data 
points  so  the  transformed  data  can  be  compressed  at  a  higher  compression  ratio. 
Dependent  on  lossy  or  lossless  criteria,  compression  can  be  implemented  using 
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entropy  coding  or  other  technique.  The  research  can  be  tasked  in  two  aspects: 
transform  and  compression  techniques. 

Wavelet  transform  has  recently  been  found  to  be  a  powerful  tool  in  signal 
processing,  and  in  multi-resolution  analysis  in  particular.  For  multi-resolution 
applications,  wavelet  transform  filters  the  data  into  a  set  of  bands  that  correspond 
to  different  levels  of  details.  The  transformed  data  are  then  transmitted  in  bands 
from  lower  to  higher  resolution,  and  a  perfect  reconstruction  of  the  original  data 
can  be  achieved  with  all  the  data  bands.  The  issues  involved  include  selection  of 
good  wavelet  filters,  optimization  of  block  partition,  and  integer/reversible 
transform. 

For  multi-resolution  signal  processing,  wavelet  transform  serves  as  a  digital 
filter  with  compact  support.  The  aim  here  is  to  have  a  filter  that  can  maximally 
pack  the  energy  in  a  minimal  number  of  coefficients,  and  be  computationally  less 
expensive.  The  efficient  use  of  wavelet  transforms  requires  that  the  dimensions  of 
the  data  are  the  powers  of  certain  integers.  For  fast  data  manipulation,  the  data 
dimension  is  required  to  be  dyatic,  i.e.,  power  of  2.  Unfortunately,  in  the  real 
world,  the  data  dimensions  are  quite  arbitrary.  Therefore  data  partition  is  often 
needed  to  divide  the  data  into  small  blocks,  e.g.,  4x4x4  or  8x8x8.  Data 
partition  helps  eliminate/reduce  the  overhead  introduced  by  padding  when 
treating  the  data  volume  as  a  single  block.  However,  it  also  affects  the  efficiency 
of  wavelet  transform.  For  example,  if  the  data  is  smooth,  a  large  block  size  for 
wavelet  transform  will  result  in  larger  energy  packed  in  a  smaller  number  of 
coefficients.  Therefore,  to  achieve  good  wavelet  transform  effects,  data  partition 
must  be  fine  tuned. 

The  integer  wavelet  transform  technique  deals  with  lossless  compression. 
General  wavelet  transform  will  create  transform  coefficients  as  floating  numbers, 
no  matter  whether  the  original  data  are  integer  or  floating.  Lossless  compression 
of  floating  number  is  difficult.  In  practice,  a  floating  number  is  often  separated 
by  its  exponent  and  mantissa  and  compressed  separately.  Since  wavelet  transform 
decorrelates  the  data  to  increase  the  signal  entropy,  this  decorrelation  may  not  be 
carried  into  the  exponent  and  mantissa.  Using  integer  wavelet  transform  on  the 
exponent  and  mantissa  integers  of  the  original  data  may  solve  the  problem  at  the 
cost  of  double  wavelet  transform  computations.  This  extra  computational  load 
can  be  reduced  by  the  integer  operation  of  the  integer  wavelet  transform. 

Efficient  compressing  wavelet  coefficients  requires  investigation  of  the 
entropy  distributions.  For  entropy  coding,  multiple  codebooks  are  suggested  for 
different  bands  of  coefficients.  In  addition,  embedded  coding  techniques,  which 
code  bits  instead  of  bytes,  will  be  studied  to  explore  multi-resolution. 
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OTHER  RESEARCH  INTEREST:  TEXTURE  MAPPING 
FOR  FLOW  VISUALIZATION 
by  Dr.  Zhifan  Zhu 


OBJECTIVE: 

Visualizing  flow  fields  is  a  challenging  task.  Effective  and  realistic 
visualization  requires  a  vivid  display  of  both  flow  directions  and  amplitudes. 
Texture  mapping  is  a  technique  used  in  computer  graphics  to  map  a  texture 
pattern,  often  an  image,  to  an  object  surface.  The  use  of  it  to  visualize  numerical 
flow  fields  have  been  proposed.  The  idea  is  to  create  texture  images  that  map  the 
flow  fields,  and  then  animate  the  texture  images.  The  mapped  texture  image 
contains  pixel  clusters  (randomly  positioned  in  the  background)  that  are  blurred 
and  deformed  in  the  flow  directions,  and  are  shaped  by  the  flow  amplitude  to 
achieve  a  realistic  visual  effect.  Such  an  image  is  repeatedly  computed  with  a 
certain  phase  shift  applied  to  the  pixels  in  the  flow  directions  to  create  a  sequence 
of  texture  images  for  animation. 

Two  popular  texture  mapping  techniques  are  Line  Integral  Convolution 
(LIC)  and  Spot  Noise.  The  LIC  method  computes  local  streamlines  for  each  pixel 
position  of  a  random  input  image  and  determines  the  output  pixel  value  by 
integration  along  the  streamline.  The  Spot  technique  takes  an  input  image 
consisting  of  random  spots  in  the  background  and  computes  a  local  stream- 
surface  for  each  spot,  and  then  deforms  the  spots  according  to  the  local  stream- 
surface.  Since  the  Spot  method  uses  2-D  deformation,  flow  amplitudes  can  be 
visualized  via  the  varying  width  of  the  spot,  but  generally  it  is  more 
computationally  expensive  than  LIC.  Both  the  LIC  and  the  Spot  Jilgorithms  have 
been  used  in  computational  fluid  dynamics  (CFD)  applications  to  visualize  surface 
flows  near  the  objects  in  CFD  fields. 

APPROACH: 

A  frequency-modulation  (FM)  technique  was  used  in  computing  the  texture 
patterns  modulated  by  the  underlying  flow  field.  The  algorithm  is  mathematically 
straightforward  and  easily  implemented.  Similar  to  LIC,  the  FM  approach  is  a 
1-D  operation  along  the  streamline  direction;  however,  it  uses  global  instead  of 
local  streamlines,  and  convolution  is  replaced  by  frequency  modulation.  The 
technique  is  omputationally  efficient  and  produces  effective  visualization.  The  FM 
signal  can  be  mathematically  written  as: 

y(t)  =  A  cos[  ^(r) ]  =  A  cos  [cot  +  X9( t)  +  ^] 
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where  A,  (O ,  C,  and  P  are  the  time-invariant  parameters  of  the  signal 
representing  the  amplitude,  carrier  frequency,  modulation  index,  and  a  constant 
phase  shift.  The  instantaneous  frequency  of  the  FM  signal  is  the  derivative  of  the 
angle,  y\r(t): 


d\|/(t)  djtot  +  Ce(t)  +  p  d9(t) 

(Jt  dt  dt 

The  frequency,  co,  is  being  modulated  by  the  time-varying  signal  x(t),  given  by; 
x(t)  =  or  0(t)  =  j  x(s)ds. 

The  modulation  index  governs  the  magnitude  of  the  modulating  term  which, 
in  most  applications,  is  smaller  than  the  carrier  frequency.  Such  an  FM  signal 
will  have  its  instantaneous  frequency  varying  about  its  carrier  frequency.  A  large 
modulation  index  can  cause  the  signal  frequency  to  deviate  from  its  carrier 
frequency  to  an  extent  that  it  may  distort  the  signal.  This  is  known  as  over¬ 
modulation. 

The  effect  created  by  texture  mapping  the  flow  field  requires  that  input 
textures  are  deformed  along  the  flow  directions  such  that  the  textures  become 
enlongated  with  larger  flow  amplitudes  and  shortened  with  smaller  amplitudes. 
This  effect  can  be  achieved  by  using  an  FM  signal  to  map  the  textures  as  follows. 
Note  that  the  value  of  an  FM  signal  alternates  periodically,  and  the  rate  of  change 
is  dependent  on  its  instantaneous  frequency.  Thus,  if  the  frequency  is  modulated 
by  the  flow  amplitude,  the  amplitude  of  the  FM  signal  will  exhibit  the  same 
variation  pattern  as  that  of  the  flow  speed.  It  is  desirable  that  the  larger  speeds  be 
mapped  into  enlongated  textures.  In  terms  of  the  FM  signal,  this  corresponds  to  a 
slower  rate  of  change.  Thus,  we  need  to  modulate  the  FM  frequency  with  a  time- 
varying  function  which  is  inversely  proportional  to  the  speed. 

The  algorithm  can  be  described  by  the  following  pseudo-code: 

Inputs:  N  streamlines  of  various  length  computed  in  the  image  domain  with 

instantaneous  speed  (flow  amplitude)  recorded  at  each  node,  carrier 

frequency  w,  and  modulation  index  C; 

Normalize  the  speed; 

for  j=l  to  N  do 

Compute  colors  for  each  node  k  as 

colorlj][k]  =  cos(wt  +  C(l-speed[j][k])  +  b).0.5-i-l. 
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end  for 

where  t  is  the  line  integral  from  the  node[j][0]  to  node[j][k]  and  p  is  a 
random  initial  phase  in  [0,27r]. 

The  texture  color  is  a  mapping  of  the  FM  signal  amplitude  and  varies 
between  0  and  1.  Since  the  algorithm  computes  the  texture  for  each  pixel  along 
streamlines  according  to  the  FM  signal  without  image  composition,  the  input 
streamlines  are  required  to  cover  the  whole  image  domain  without  overlapping. 
Otherwise,  signal  interference  will  occur  and  damage  the  patterns.  On  the  other 
hand,  to  obtain  the  best  modulation,  the  streamlines  should  be  as  long  as  possible. 
To  meet  this,  the  following  algorithm  is  used  to  generate  the  streamlines: 

Inputs:  flow  field,  minimum  length  L 

Sub-divide  the  flow  field  into  small  cells  such  that  each  cell  corresponds  to  a 
pixel  in  the  output  image; 

Create  streamlines  starting  from  each  cell; 

Sort  the  streamlines  from  maximum  length  to  minimum  length  and  delete 
streamlines  that  shorter  than  L  ->  N  streamlines  of  at  least  L  long 
Reset  cell  mask  to  zeros 
for  j=l  to  N  do 

Check  the  cells  that  the  streamlinelj]  passes 

if  no  cells  marked,  then  mark  these  cells  traveled  by  streamlinelj] 

else  if  all  cells  already  marked,  discard  streamline[j] 

else  break  streamlinelj]  into  small  streamlines  according  to  the  cell  mask, 
discard  those  shorter  than  L  and  insert  the  rest  into  the  streamline  array 
according  to  their  length,  delete  streamlinelj] 
and  adjust  N 
end  for 

An  adaptive  time  increment  ‘nf  is  used  in  the  streamline  integration  to 
ensure  the  distance  between  adjacent  nodes  of  the  streamline  is  about  the  same 
length  and  cell  size.  Thus,  the  streamline  length  can  be  approximated  by  the 
number  of  nodes.  The  final  cell  coverage  depends  on  the  minimum  length  L.  This 
constraint  is  used  to  eliminate  short  streamlines  that  will  not  have  the  contribution 
to  the  final  texture.  If  L  is  chosen  to  be  one,  all  cells  will  be  covered. 

RESULTS: 

Figure  8  shows  results  from  three  different  flow  fields  depicted  at  different 
resolutions.  Figure  8a  depicts  a  synthetic  circular  flow  field  of  280  x  280  pixels. 
The  flow  amplitude  monotonically  increases  from  the  center  to  the  outside. 
Figure  8b  shows  a  horizontal  wind  field.  The  original  flow  field  has  a  lower 
resolution  than  the  image  (264  x  210).  The  larger  wind  speed  is  in  the  middle  of 
the  field  causing  the  enlongated  texture  patterns.  Two  cyclone  centers  are  easily 
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visible.  Figure  8c  illustrates  the  top  layer  output  of  an  ocean  circulation  model  of 
the  Gulf  of  Mexico.  It  features  large  ring  structures  and  many  small  turbulent 
flows.  Since  it  is  merely  a  demonstration  of  the  effectiveness  of  our  FM-based 
algorithm,  we  chose  not  to  perform  any  subdivision.  Each  pixel  corresponds  to  a 
cell  in  the  field  with  194  x  146  cells.  A  minimum  length  of  5  was  used  in  creating 
the  streamlines  for  the  wind  and  the  current  fields,  resulting  in  a  small  number  of 
pixels  not  being  covered  by  the  streamlines,  but  instead  were  assigned  the 
background  color. 
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Figure  8.  FM  Texture  mapping:  a)  a  synthetic  circular  flow  field  of  280  x  280, 
b)  a  horizontal  wind  field  of  264  x  210,  and  c)  an  ocean  circulation  field  of 
194  X  146. 


41 


Comparing  the  results,  the  circular  flow  field  has  the  best  effect.  This  is 
because  the  field  is  a  perfect  circular  flow.  There  is  no  flow  convergence  or 
divergence.  So,  except  for  the  four  comers,  all  streamlines  are  complete  circles. 
The  texture  of  the  wind  field  also  displays  a  quality  flow  image.  The  textures  in 
the  middle  exhibit  enlongated  patterns  due  to  the  larger  wind  velocities;  the  areas 
with  slower  speeds  at  the  upper-left  comer  and  near  the  bottom  have  shorter 
texture  patterns.  However,  some  discontinuities  can  be  found  in  the  areas  with 
convergence  and  divergence.  Particularly  in  regions  where  divergence  occurs  (to 
the  right  of  the  upper  cyclone  center),  such  discontinuities  are  easily  noticeable. 
This  is  because  many  new  streamlines  are  started  there  to  fill  the  decompressed 
area,  and  the  texture  colors  at  the  beginning  of  those  streamlines  are  decided  by 
the  random  phases,  producing  abmpt  effects.  The  texture  image  for  the  ocean 
circulation  field  has  the  lowest  resolution  in  the  three  examples  but  much  more 
flow  turbulence.  As  a  result,  the  created  streamlines  consist  of  more  short  pieces, 
and  the  minimum  length  constraint  left  relatively  larger  portion  of  pixels 
uncovered.  Even  so,  the  texture  is  still  able  to  show  the  flow  patterns  of  the 
eddies  with  various  flow  amplitudes. 

FUTURE  RESEARCH: 

The  main  computation  of  the  algorithm  is  in  sorting  and  reorganizing  the 
streamlines  to  eliminate  overlap  and  to  keep  them  with  maximum  possible 
lengths.  The  current  algorithm  used  to  create  such  streamlines  is  virtually  a 
brute-force  method,  and  further  efforts  are  required.  Once  the  global 
streamlines  are  generated,  the  mapping  is  quick.  This  allows  easy  modification  of 
mapping  parameters,  and  hence  gives  the  algorithm  the  desired  flexibility  and 
control.  Further  research  will  also  be  directed  at  how  to  smooth  the  signal 
abruptness  at  the  two  ends  of  the  streamlines.  Such  abruptness  occurs  because  the 
signal/texture  is  only  continuous  within  the  streamline,  and  the  signal  at  the  first 
point  of  a  streamline  is  determined  by  a  random  phase  b.  This  is  particularly 
problematic  in  large  compression  and  decompression  areas. 
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